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Abstract. 


Samples  of  biological  tissue  are  modelled  as  inhomogeneous  fluids 
with  density  p(x)  and  sound  speed  c(x)  at  point  x.  The  samples  are 
contained  in  the  sphere  |x|  <6  and  it  is  assumed  that  p(x)  =  p0  =  1  and 
c(x)  =  c0  =  1  for  |x|  >  5,  and  (Yn(x) j  <<  1,  [Yp(x) j  «  1  and 
|Vyp(x) j  «  1  where  yp(x)  =  p(x)  -  1  and  Yn(x)  =  c  (x)  -  1.  The  samples 
are  insonified  by  plane  pulses  s(x  •  0O  -  t)  where  |0O|  =  1  and  the 
scattered  pulse  is  shown  to  have  the  form  jx|  1  es(|x|  -  t,0,0o)  in  the 
far  field,  where  x  »  |x|0.  The  response  es(T,0,0Q)  is  measurable.  The 
goal  of  the  work  is  to  construct  the  sample  parameters  yR  and  yp  from 
es(T,0,0o)  for  suitable  choices  of  s,  0  and  0O • 

In  the  limiting  case  of  constant  density:  Yp(x)  =  0  it  is  shown 

that 


Vx>  -  i 


:^(2x  •  9, 0,-0)  d0 


where  5  represents  the  Dirac  5  and  S2  is  the  unit  sphere  | © [  =  1. 
Analogous  formulas,  based  on  two  sets  of  measurements,  are  derived  for 
the  case  of  variable  c(x)  and  p(x) . 


ii 


1.  Introduction . 


Computer-based  acoustic  imaging  techniques  have  been  studied 
intensively  during  the  last  decade  [6 ,8,9]*  Typical  techniques  involve 
irradiating  a  sample  with  prescribed  sound  fields,  measuring  the 
resulting  scattered  fields  and  applying  a  computational  algorithm  to  the 
scattering  data  to  produce  maps  of  such  sample  parameters  as  density, 
sound  speeds  and  perhaps  others.  These  techniques  have  important  appli¬ 
cations  to  medical  ultrasonic  imaging  where  the  sample  is  a  living 
organism,  to  non-destructive  evaluation  where  the  sample  is  a  manufactured 
item  such  as  a  metal  casting  or  ceramic  object  and  to  geophysical 
prospecting  where  the  sample  is  a  portion  of  the  earth's  crust. 

This  paper  treats  a  problem  of  medical  ultrasonic  imaging.  The 
sample  is  modelled  as  an  inhomogeneous  fluid  which  is  characterized  by  a 
variable  density  p(x)  and  sound  speed  c(x) .  The  use  of  a  fluid  model  is 
motivated  by  the  fact  that  in  biological  tissues,  other  than  bone,  acoustic 
shear  waves  are  not  observec . 

In  acoustic  imaging,  the  scale  of  the  smallest  structures  that 
can  be  resolved  is  of  the  order  of  the  smallest  wavelength  employed.  A 
typical  sound  speed  in  biological  tissue  is  c  =  1500  m/sec.  Thus  for  a 
wavelength  of  A  *  1  mm.  =10  3  m. ,  a  frequency  of  f  *  c/i  «  1.5  *  106  hertz 
=  1.5  megahertz  is  needed.  This  is  in  the  high  ultrasonic  range. 

The  acoustic  field  in  an  inhomogeneous  fluid  with  density  p(x) 
and  sound  speed  c(x)  may  be  characterized  by  a  real  valued  function 

♦Numbers  in  square  brackets  denote  references  from  the  list  at  the  end 
of  the  paper. 
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u(t,x)  that  satisfies  the  scalar  wave  equation 

(1.1)  Ijpr  -  C2(x)  p(x)  V  •  7u  *  0 

where  t  is  a  time  coordinate  and  x  =  (xj.X2.X3)  €  R3  denotes  a  triple  of 
rectangular  coordinates  in  space.  u(t,x)  may  be  interpreted  as  the 
acoustic  pressure;  that  is,  the  difference  between  the  instantaneous 
pressure  and  the  equilibrium  pressure  in  the  fluid.  It  is  directly 
measurable.  Derivations  of  (1.1)  from  the  principles  of  fluid  dynamics 
may  be  found  in  a  number  of  books  and  monographs;  see,  e.g.  [1,2, 3, 7]. 

It  will  be  assumed  here  that  the  sample  to  be  imaged  is 
contained  in  the  ball  B(0,6)  =  (x  :  )x|  <  6},  and  that 

(1.2)  p(x)  5  pa  =  1  and  c(x)  =  c„  »  1  outside  B(0,6), 

where  p0  and  c0  are  the  constant  parameters  of  the  ambient  fluid.  The 
conditions  P0  *  1,  c0  *  1  are  a  convenient  normalization  that  can  always 
be  achieved  by  a  suitable  choice  of  units. 

The  deviations  of  the  sample  parameters  from  those  of  the  ambient 
fluid  will  be  measured  by  the  parameters 

(1.3)  Yp(x)  -  p(x)  -  1 
and 

(1.4)  Yn(x)  ”  c“2(x)  -  1  =  nz(x)  -  1 

where  n(x)  ■  c  1 (x)  is  the  index  of  refraction.  The  acoustic  imaging 


method  developed  below  is  based  on  the  Bim  approximation  to  solutions 
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of  (1.1).  The  conditions  for  the  validity  of  the  approximation  are 

(1.5)  |yn(x) I  «  1,  !yp(x) i  «  1  and  |Vyp(x)|  «  1. 

The  purpose  of  an  ultrasonic  imaging  technique  for  (1.1)  is  to 
construct  accurate  maps  of  the  functions  p(x)  and  c(x)  by  applying  a 
computational  algorithm  to  the  measured  data  of  suitable  scattering 
experiments.  An  explicit  method  for  solving  this  two  parameter  imaging 
problem  is  derived  below  under  the  weak  scattering  hypothesis  (1.5). 

The  basic  scattering  experiment  is  the  scattering  by  the  sample  of  a 
plane  wave 

(1-6)  u0(t,x)  =  s(x*  60  -  t) 

where  9o  is  a  fixed  unit  vector,  or  point  on  the  unit  sphere 
S2  *  {x  :  | x |  =  l)  C  R3,  and  s(t)  is  a  prescribed  wave  profile. 

The  imaging  method  developed  below  is  based  on  the  author's 
theory  of  asymptotic  wave  functions  as  developed  in  [14,15,16].  In  the 
present  context  the  theory  states  that  if  u(t,x)  is  the  total  field  due 
to  the  interaction  of  the  pulse  (1.6)  with  the  sample,  and  if 

(1.7)  uSc(t,x)  =  u(t,x)  -  u0(t,x) 

sc 

is  the  scattered  field  then  u  has  the  far  field  form 

(1-8)  uSC(t,x)  «  | x |  1  eg(|x|  -  t,9,90)  +  o(l) , 

where  x  ■  |x|8,  9  €  S2  and  the  error  term  o(l)  tends  to  zero  when  t 
The  imaging  algorithm  developed  below  is  based  on  the  echo 
profile  function  e  (t,6,0o).  The  results  take  their  simplest  form  when 
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s(t)  «  5(t),  the  Dirac  6-function.  Of  course,  an  ideal  infinitely  sharp 
pulse  cannot  be  realized  in  practice.  However,  good  approximations  can 
be  generated.  Alternatively,  one  can  obtain  the  response  6^(1, 9, 90)  by 
electronic  filtering  of  the  response  es(T,0,0o)  to  a  more  realistic  pulse 
profile  s(t). 

The  imaging  method  takes  a  particularly  simple  form  in  the 
constant  density  case  (yp(x)  =  Q)  when  there  is  only  one  parameter  to  be 
imaged.  In  this  case  it  will  be  shown  that,  in  the  Born  approximation. 


(1.9) 


Yn(x) 


e^(2x  •  0,9, -0) 


d9 


where  d0  is  the  element  of  area  (solid  angle)  on  S2.  Thus  is  obtained 
by  integrating  the  back  scattered  echoes  eg(T,0,-0)  over  all  directions. 
In  the  general  case  a  second  set  of  measurements  is  needed  to  determine 
the  parameters  Yp  and  Yn- 

A  formula  equivalent  to  (1.9)  was  derived  by  S.  J.  Norton  and 
M.  Linzer  [9]  who  obtained  it  as  a  limiting  case  of  an  imaging  method 
based  on  near  field  measurements  (see  [9,  p.  215,  (81)]).  More  recently, 
J.  H.  Rose  and  J.  M.  Richardson  [11]  have  given  without  proof  an  analogue 
of  (1.9)  for  the  imaging  of  inhoraogeneities  in  isotropic  elastic  solids. 
In  their  work  they  also  formulate  analogues  of  (1.8)  for  elastic  solids 
and  discuss  their  applications  to  multiparameter  imaging. 

The  two  goals  of  this  paper  can  now  be  formulated.  The  first 
goal  is  to  show  how  the  author’s  theory  of  asymptotic  wave  functions, 
cited  above,  and  the  Born  approximation  can  be  used  to  provide  a  simple 
and  rigorous  derivation  of  a  functional  relation  between  the  echo  wave¬ 
form  e^x, 0,0O)  and  the  Radon  transforms  of  Yn  and  Yp-  The  second  goal 
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is  to  use  this  relation  to  derive  an  explicit  imaging  algorithm  for  the 
two  parameter  case.  It  will  be  seen  from  the  analysis  that  the  same 
method  is  applicable  to  electromagnetic  imaging,  and  to  acoustic 
imaging  of  solids  where  more  than  two  parameters  are  to  be  imaged. 

The  remainder  of  the  paper  is  organized  as  follows.  §2  presents 
a  brief  discussion  of  the  facts  concerning  the  scattering  of  time- 
harmonic  plane  waves  that  are  needed  to  analyze  pulse  mode  scattering. 
This  theory  is  similar  to,  and  simpler  than,  the  theory  of  scattering  by 
bounded  obstacles  developed  in  [14,16].  §3  presents  the  Born  approxima¬ 

tion  to  the  time-harmonic  scattered  fields  and  the  scattering  amplitude. 
§4  develops  the  functional  relation  between  the  scattering  amplitude  and 
the  Radon  transforms  of  yn  and  y  .  §5  reviews  the  theory  of  pulse  mode 

scattering  as  a  boundary  value  problem  and  the  associated  theory  of 
asymptotic  wave  functions.  The  final  §6  develops  the  acoustic  imaging 
method  for  the  two  parameter  problem.  For  clarity  the  known  one 
parameter  formula  (1.9)  is  derived  first.  Then  the  method  is  extended 
to  the  general  two  parameter  case.  The  section  ends  with  a  brief 
discussion  of  the  numerical  implementation  of  the  method. 


2 .  The  Scattering  of  Time-Harmonic  Plane  Waves. 

A  time-harmonic  plane  wave  propagating  in  the  direction  0Q  E  S2 
in  the  ambient  fluid  may  be  characterized  by  the  complex  wave  function 

(2.1)  w #(x,u)0o)  =  <2tt) “3/2  eiw0°‘x 

where  gV2it  is  the  wave  frequency  and  the  time  factor  e  has  been 
suppressed.  The  amplitude  factor  (2 n)  3^2  is  a  normalization  that  is 
included  to  facilitate  the  application  below  of  the  results  of  [16]. 

If  the  scatterer  is  irradiated  by  the  field  (2.1)  the  resulting 
time-harmonic  field 

(2.2)  w(x,w60)  =  w0  (x,u)0o)  +  wSc(x,ue0) 
is  uniquely  characterized  by  the  properties 


(2.3) 


c2 (x)  p(x)  V  • 


1 

P(x) 


Vw 


w2w 


0 


for  all  x  G  R3  and 


(2.4) 


3wsc 
3  { x  | 


iujwsc  =  0(|x|  2),  |xj  -*■<*>  . 


Equation  (2.3)  is  just  the  wave  equation  (1.1)  for  w(x,a)90)  e  iwt.  (2.4) 
is  the  Sommerfeld  outgoing  radiation  condition.  (2.3)  and  (2.4)  will  be 

SC 

shown  to  imply  that  w  has  the  far  field  form 


(2.5) 


wsc(x,w0o) 


iujxj 

"4if[xJ  T(w8>w9o)  +tf(M~2),  M  -  06 , 
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where  x  *  | x } 8 .  The  coefficient  T (p , p ' )  is  called  the  scattering 
amplitude  of  the  scatterer.  It  is  known  to  satisfy  the  reciprocity  law 


(2.6)  T(p.p’)  =  T(-p\-p). 

T(p,p')  plays  a  key  role  in  the  theory  of  acoustic  imaging  presented 
below. 

The  wave  field  w(x,p)  and  its  scattering  amplitude  can  be 
constructed  by  solving  an  integral  equation.  To  derive  it  note  that 

(2.7)  V  •  wj  =  w  -  Vp  •  7w  . 

Thus  (2.3)  may  be  written 

2 

(2.8)  Aw  -  (V  £n  p(x))  •  7w  +  g'Zy^y  w  «  0  . 

This  can  be  treated  by  perturbation  theory  by  introducing  the  parameters 

Y  and  Y. . 
n  P 

Note  that  by  (1.3) 

(2.9)  V  in  p  =  V  J,n  (1  +  yp)  *  (1  +  Yp)_1  V  Yp- 
Hence  (2.8)  is  equivalent  to 

(2.10)  Aw  +  oj2w  ■  -w2  Yn(*)w  +  (1  +  Yp(x))-1  Vyp(x)  •  Vw 
by  (1.4)  and  (2.9).  Recall  that  by  hypothesis 

(2.11)  F  =  supp  Yn  U  supp  Yp  C  B(0 ,5) 

where  supp  Yn  denotes  the  support  of  Yn  (the  smallest  closed  set  outside 


of  which  Y  (x)  =  0).  Application  of  Green's  theorem  to  w(x,p)  and  the 
outgoing  Green's  function  for  the  ambient  fluid 


(2.12) 


G(x-x',p)  = 


l|p||X-X 


4tt  I  x-x' 


gives,  after  a  standard  calculation, 
w(x,p)  =  w0(x,p) 

(2.13) 


+ 


G(x-x' , p) { | p | 2Yn(x' )  w(x' ,p)  -  (1+Yp(x') ) 


Vy.(x')  •  Vw(x' ,p)}  dx' 

K 


where  dx'  =  dxj  dx^  dxj .  If  x  is  restricted  to  T  then  (2.13)  is  an 
integro-dif f erential  equation  for  w(x,p)|F,  the  restriction  of  w(x,p)  to 
F.  It  may  be  transformed  into  a  pure  integral  equation  by  an  integration 
by  parts  in  the  last  term.  Solution  of  this  integral  equation  by 
standard  techniques  provides  a  construction  of  w(x,p)|r.  The  continua¬ 
tion  of  w(x,p)|r  to  all  of  R3  is  then  provided  by  (2.13). 

A  verification  of  (2.5)  and  a  construction  of  T(p,p')  are  also 
provided  by  (2.13).  Indeed,  if  x  =  ] x | 0  then 

(2.14)  |x-x'|  =  j  x  |  -  x'  •  6  +  0(  |x|  *)  ,  !  x  |  -*■  », 

uniformly  for  all  x'  €  f.  Thus 


(2.15)  G(x  -  x\p) 

uniformly  for  x’  6  f. 


-i  j  p  j  0 -x  ’ 

4tt  |x  I  6 


Substitution  of  (2.15) 


+  0( jx [  2) , 

into  (2.13) 


gives  (2.5) 


with 


10 


T(u)0,oj0o)»  e-iw0*x  {^^(x)  w(x,u>0o)  -  (l-^D(x))-17Yp(x)  *7w(x,u)0o) }  dx 

(2.16) 

Thus  T  can  be  calculated  from  the  parameters  Yn>  Yp  and  the  field  w(x,p) 
inside  the  scatterer. 


3.  Weak  Scatterers  and  the  Born  Approximation. 

When  the  scatterers  are  weak  in  the  sense  of  conditions  (1.5)  the 
integro-dif f erential  equation  (2.13)  for  w(x,p)  can  be  solved  by 
iteration.  On  dropping  terns  of  orders  higher  than  the  first  in  Yn»  Yp 

SC 

and  Vyp  one  obtains  the  Born  approximation  to  w  (x,p).  It  is  given  by 


sc ,  „ 

w  (x,p) 

(3.1) 


G(x-x',p)  { 1 P  |  2  Yn (x ' )  w0(x',p)  -  Vy0(x’)  •  Vw0 (x ' ,p) )  dx ' 


(2ir)  1/2  I  G(x-x',p)  eip  x 


{|p!2Yn(x’)  -  ip  •  VYp(x')}  dx’. 


The  corresponding  Born  approximation  to  the  scattering  amplitude  is 


(3.2)  T(a)6,we0)  =  (2tt) 


r 

n 

■  r 


3/2  [  o  iw(.-0o)  X  {u2-y  (x)  -  iw6a  *  Vy.(x)}  dx. 


An  integration  by  parts  in  the  last  term  gives  the  alternative 
representation 


(3.3)  T(m8.u)90)  -  j  e'iw(9"0o) ’x  (Yn(x)  +  (0  •  90  -  1)  Yp(x) }  dx. 


The  notation 

(3.4)  Yp(x)  “  Yfl(x)  +  (y-1)  Yp(x) 

will  be  used  in  what  follows.  Equation  (3.3)  can  then  be  written 
concisely  as 
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(3.5)  T (w9,u)90 )  -  w2  Yg.g^  (w(9  -  80)) 


where  y  is  the  usual  Fourier  transform,  defined  by 


(3.6) 


Y(P)  =  (2-rr) 


■3/2 


e  ip  x  y(x)  dx. 


In  the  remainder  of  this  paper  the  approximation  (3.5)  will  be  used  for 
the  scattering  amplitude. 


4 .  The  Born  Approximation  and  the  Radon  Transform. 


The  Radon  transform  of  a  function  y(x)  with  compact  support  is 
the  function  Y  :  R  *  S2  ■*  R  defined  by  [5,10] 


(4.1) 


Y(s,n) 


Y(x)  dS 

x*n=s 


where  (s,n)  6  R  *  S2.  Equation  (4.1)  means  that  y(s,n)  is  the  integral 
of  Y(x)  over  the  plane  with  equation  x  •  0  =  s  and  surface  element  ds. 
An  alternative  notation  is 


(4.2) 


Y(s,n)  = 


Y(sn  +  x  )  dx' 


where  is  a  point  in  the  plane  through  the  origin  that  has  normal  n  and 

surface  element  dx-*-.  Note  that  y(s,n)  also  has  compact  support: 

supp  Y  C  B(0,5)  =*  Y(s,n)  =  0  for  |s|  >  5. 

In  this  section  a  known  relation  between  the  Fourier  and  Radon 

transforms  is  derived  and  used  to  relate  the  Born  approximation  (3.5) 

to  the  Radon  transform  of  ya  a  .  This  relation  and  the  known  inverse 

“  *  “o 

Radon  transform  provide  the  basis  for  the  imaging  method  developed  in  §6. 
To  begin,  consider  the  Fourier  integral 


(4.3)  (2tt)3/2  y(u(0- 60))  =  e_iu)(9"eo)’x  y(x)  dx 

Jr3 


where  u>,  9  and  90  are  fixed,  tu  ^  0  and  9  #  9q  •  Introduce  new  variables 
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(4.4) 


n  .  _  c  c2 

n  s  ’ 


(4.5) 


t  -  (6  -  e0)  •  x, 


and 


(4.6) 


s  =*  x  •  n 


19-00 


Then 


(4.7) 


x  =  (x  •  n)n  +  x  =■  sn  +  x  , 


where  x  is  in  the  plane  orthogonal  to  n,  and  a  rotation  of  coordinates 
in  (4.3)  gives 


(2tt)  3/2  y(uj(9  -  e0})  «  [  e~iaJT  y(srt  +  x-1-)  ds  dx~ 

JR3 


(4.8) 

=  f  e'iuT 

l  i 

Y(sri  +  x  )  dx 

J  —CO 

R2 

ds 


0  ~  0r 


1  f  e 

-iwx  - 

f  T 

9-0o  } 

L 

Y 

ll^eTT 

’  TFeTlJ 

dr. 


The  last  integral  is  a  one  dimensional  Fourier  transform.  Hence 
multiplying  it  by  u)2  is  equivalent  to  the  operation  -32/9t2  :>n  the 
integrand.  Thus 


(4.9)  (2tt)3/2  o,2y(u)(0-0o))  - -|0-0Df3  [ 

J  — OC 


e-iWT  r 


0-0o 


9-00  ’  0-9 


o  IJ 


dx 


where  Y"(s,n)  denotes  the  second  s  derivative  of  y(s»h).  Now  the  function 
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(4.10) 


h  (s,n)  =  - 


Y"(s,n) 


occurs  in  Che  inverse  Radon  transform.  In  fact,  the  inverse  is  given 
by  [5] 

(4.11)  Y(x)  =  [  h  (x  •  n,n)  dn. 

JS2  Y 

Thus  it  is  natural  to  rewrite  (4.9)  as 


(4.12)  u2  Y(u»(9-9t))  =  f-effp- 


-OO 

A*iU)T  h 

f,  T  0-00  } 

N 

-oo  1 

[|0-0O  ’  I e-0o  J 

dT . 


Applying  (4.12)  to  the  Born  approximation  (3.5)  gives  the  relation 

1/2  f°°  !  T 


(4.13)  T(co0  ,w0o)  =  t 


0-0o  |  3 


—  1UT  , 
e  h 


0*6n 


0-00 


e-eB 


0-0  „ 


dx 


where,  for  brevity,  the  notation 

(4.i4)  hQ.eo  =  h  =  Yn  +  (0  •  e0  -  l)Yp 

has  been  introduced. 

The  scattering  amplitude  can  be  obtained  from  far  field 
measurements;  see  (2.5).  Hence,  the  relation  (4.13)  is  a  natural 
starting  point  for  an  imaging  method.  The  Fourier  transform  of  (4.13)  is 


(4.15) 


1 

(2tt)  ‘/z 


r°° 

eiTU)  T(w0,w0o)  du 

—00 


4 11  h 

r 

T 

0-00 

i e-e0 1 3  e-0o 

9-00 

’  J  0—0  0 

Suitable  choices  of  r,  9  and  90  in  this  relation,  with  0  •  0O 


U  fixed. 
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give  hy(s,n) .  Then  the  inversion  formula  (4.11)  gives  "y  .  By  doing 
this  for  two  values  of  y  one  can  deduce  y ^  and  from  (3.4).  This 
program  is  carried  out  in  §6  below.  Moreover,  it  is  shown  that  the 
values  of  the  Fourier  integral  on  the  left-hand  side  of  (4.15)  can  be 
obtained  directly  from  pulse  mode  scattering  measurements.  This  is 


developed  in  the  next  section. 


5.  Pulse  Mode  Scattering. 

The  problem  of  the  scattering  by  the  sample  of  plane  wave  pulses 

(5.1)  u0(t,x)  ■  s(x  •  0#  -  t) 

Is  formulated  and  solved  In  this  section.  The  pulse  profile  s(t)  1* 
assumed  to  be  a  prescribed  function  with  compact  support: 


(5.2)  supp  s  C  [a,b]  . 

Thus  for  any  fixed  value  of  t  one  has 

(5.3)  supp  u0(t,*)  c{x:a+t<x*60<b+t) 

and  the  pulse  (5.1)  does  not  interact  with  the  sample,  which  is 

contained  B(0 ,6) ,  before  the  time  to  »  -b  -  6  :  see  Figure  1. 

Sample 


Figure  1.  Incident  pulse  ug  before  Interaction  with  the  sample. 


Therefore,  the  total  acoustic  pressure  field  u(t,x)  due  to  the  scattering 
of  the  pulse  (5.1)  satisfies 
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(5.4)  u(t,x)  =  u0(t,x)  for  all  t  <  t0  and  x  €  R3. 

It  follows  that  u(t,x)  is  a  solution  of  an  initial  value  problem  for  the 
wave  equation  (1.1): 


(5.5) 


c2(x)  p(x)  V 


=  0  for  all  t  >  t0 


and  x  6  R 3 , 


(5.6) 


u(to,x)  =  f  (x)  and  °u--  =  g(x)  for  all  x  e  R3 


where 


(5.7) 


' 

f(x) 

g(x) 


u0(t0  »x)  =  s(x  *  6o  -  t„) ,  and 
9. -t.)  . 


The  theory  of  the  initial  value  problem  (5.5),  (5.6)  is  analogous 
to,  and  simpler  than,  the  theory  of  initial-boundary  value  problems  for 
the  wave  equation,  as  developed  in  [14,15,16].  The  theory  will  be 
outlined  here  without  proofs.  Details  may  be  found  in  the  references 
cited. 

The  problem  (5.5),  (5.6)  is  most  simply  discussed  in  terms  of 
the  scattered  field 


(5.8)  usc(t,x)  =  u(t,x)  -  u0(t,x). 

It  is  a  solution  of  the  problem 


(5.9) 


*2  SC 
3  u 

"3P“ 


-  c2 (x)  p(x)  7  • 


1 

P(x) 


F(t,x)  for  t  >  t0 ,  x  s  R3 , 
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(5.10) 


sc 

u 


(t0 ,x) 


=  0  and 


5uSC (t  o ,x) 
dc 


0  for  x  €  R3 


where  F(t,x)  is  defined  by 


(5.11) 


F(t.x)  = 


82uo (t ,x) 

— aP — 


+  c2 (x)  p (x)  V 


P(x) 


Vu,(t, 


x)j 


for  all  t  >  t0  and  x  £  R3 .  Note  that  (5.3)  and  the  assumption  that 
p(x)  =  1,  c(x)  *  1  outside  B(0,5)  imply  that  F  has  compact  support  in 
space-time.  More  precisely, 

(5.12)  supp  F  c  [ -b  —  <5]  *  B(0,5)  . 

A  simple  approach  to  the  initial  value  problem  (5.9),  (5.10)  may  be 
based  on  the  theory  of  the  operator 


(5.13) 


Au 


-c2(x)  p(x)  V 


1 

P(x) 


Vu 


in  the  Hilbert  space 

(5.14)  JC  =  L2(R3,c'2(x)  p-1 (x)  dx) 
with  scalar  product 

(5.15)  (u,v)  »  f  u(x)  v(x)  c  2(x)  p  1 (x)  dx  . 

JR3 

The  theory  may  be  based  on  Kato's  theory  of  sesquilinear  forms  in  Hilbert 
spaces  [4];  see  [17] .  It  follows  that  if  the  domain  of  A  is  defined  by 
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(5.16) 


D(A)  -  3C  n 


Vu  and  V 


IVu 


are  in 


* } 


Chen  A  is  a  selfadjoint  non-negative  operator  in  The  problem  (5.9), 
(5.10)  may  then  be  reformulated  as  the  Hilbert  space  problem 


(5.17) 


2  SC 

d^Z~  +  AuSC  =  F(t ,  •)  for  t  >  t0, 


(5.18) 


uSC(t0)  =  0  and 


duSC(to) 

dt 


The  formal  solution  is  given  by  the  Duhamel  integral 


(5.19) 


uSC(t,*) 


t 1 


{A  l^2  sin  (t-t)  A1/2}  F(t,*)  dT  . 


1/2’ 


A  rigorous  interpretation  of  (5.19)  may  be  based  on  the  spectral  theorem 
for  A.  In  particular,  uSC  is  a  "strict  solution  with  finite  energy"  in 
the  sense  of  [12]  provided  that  t  -*•  F(t,')  e  3f  is  continuous.  If  p(x)  , 
c(x)  and  s(x)  are  smooth  functions  then  known  regularity  results  [13] 
imply  that  u  is  a  classical  solution. 

It  will  be  convenient  here  to  represent  usc(t,x)  as 

(5.20)  usc(t,x)  *  Re  (vsc(t,x)} 

SC 

where  v  is  the  complex  valued  wave  function  defined  by 


(5.21) 


vsc  (t ,  •)  •  i  A 


■1/2 


-itA 


1/2  ft 


ItA1/2  ,  . 

e  F(x,*)  dt  . 


Note  that,  since  F(t,*)  •  0  for  t  >  tj  =  -a  +  S,  one  has 

bc  -itA1/2 

(5.22)  vsc(t,*)  »  e  h  for  t  >  tx 
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where 


(5.23) 


h  *  i  A 


‘1/2 


ftl  i?A1/2  ,  .  , 

e  F(t,«)  dr 


The  asymptotic  behavior  for  t  -*■  00  of  wave  functions  of  the  form 
(5.22)  was  calculated  in  [14,16]  for  solutions  of  the  d'Alembert  equation 
in  exterior  domains.  By  the  same  methods  it  can  be  shown  that  if 
x  *  | x | 0  then 

(5.24)  vSC(t,x)  =  |xflF(  |x|  -  t,0)  +0(1) 


where 


(5.25) 


dQdx  <  00 


and  0(1)  -*■  0  in  5C  when  t  The  function 


(5.26) 


v®C(t,x)  =  Ixf1  F(  ]x  |  -  t,9) 


s  c 

is  called  the  asymptotic  wave  function  for  v 

The  results  (5.24),  (5.25)  are  based  on  the  fact  that  if 


(  w+(x,p)  =  w(x,p) 


(5.27) 


w_(x,p)  =  w+(x,-p)  =  w(x,-p)  J 


for  x  €  R3,  p  €  R3, 


where  w(x,p)  is  the  time-harmonic  field  of  §2,  then  each  of  the  families 
(w+(x,p)  :  p  6  R3}  and  (w_(x,p)  :  p  e  R3}  is  a  complete  family  of 
generalized  eigenfunctions  for  the  operator  A.  In  fact,  if 
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These  results  may  be  verified  by  the  method  of  [14]  and  [16]. 

By  taking  the  real  part  of  (5.24)  one  gets  the  asymptotic  form 
of  usc(t,x)  for  large  t: 

(5.32)  uSC(t,x)  -  jx|  1  eg(  |x|-t,0,0o)  +0(1) 

where  eg(T,0,0o)  ”  Re  {F(t,9)}  and  0(1)  ■*  0  in  X  when  t  ->  ®.  By  (5.30) 
the  wave  profile  eg  is  given  by 


e  (t,9,90) 


e1Ta)  T(uj8 ,u)8n)  s (to)  du> 


eiTU)  T((i)0,a)9o)  §(u)  du  . 


(5.33) 
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The  last  equation  follows  from  the  assumption  that  s(x)  is  real,  whence 
s(-x)  =  s(t) ,  and  the  analogous  property  of  the  scattering  amplitude: 

(5.34)  T  ( -0)6 ,  -u)0  0 )  =  T(u)0,a)0o)  . 

Relation  (5.32)  implies  that  the  pulse  echo  profile  es(t,0,0o) 
may  be  obtained  from  far  field  measurements: 

(5.35)  es(t,0,8o)  ~  |x)  uSC( |x)-t, |x|0) , 

with  an  error  that  tends  to  zero  when  t  or  equivalently  when 

]x|  •*■  *>.  This  result  and  the  relation  (5.33)  are  used  in  §6  to  construct 

an  imaging  method. 


6 .  Pulse  Mode  Imaging. 

In  this  section  the  Born  approximation  (A. 13)  and  the  representa¬ 
tion  (5.33)  for  the  pulse  echo  profile  are  combined  to  obtain  an  explicit 
relation  between  e  and  the  function  • 

First,  note  that  (4.13)  guarantees  that,  in  the  Born  approxima¬ 
tion,  T(u)6,(u0o)  is  the  Fourier  transform  of  the  well-behaved  function  of 
T  given  by  (4.15).  It  follows  from  (4.15),  (5.33)  and  the  convolution 
theorem 


(6.1) 


r°° 

f(T')  g ( T  -  T’)  dT’ 

_00 


(CO 

eiTU  f(to)  g(w)  du) 

— oo 


that 


(6.2) 


Gs(x,0,eo) 


2~tt 

1 0-00  1 5 


r  sc. 


—00 


Ve0 


T-T1 

.T0-0oT 


Note  that  on  taking  s(t)  =  5(t)  ,  the  Dirac  delta,  one  gets 


(6.3) 


i6(T,9,e0) 


2*  h 

T 

0_eo  ] 

I e-90 i 3  e-0o 

O 

CD 

1 

CD 

’  Te^eTTJ 

Relation  (6.3)  provides  the  basis  for  the  imaging  method  of 
this  paper.  The  pulse  echo  profile  e^  will  be  regarded  as  obtainable 
from  scattering  measurements.  Equation  (6.2)  shows  that  good  approxi¬ 
mations  to  e^  can  be  obtained  with  incident  pulses  s(t)  that  are  smooth 
approximations  to  6(t) .  Moreover,  as  mentioned  in  the  introduction,  e^ 
can  be  obtained  by  suitable  filtering  of  the  echoes  from  other 
profiles  s(t) . 
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The  imaging  of  the  parameters  Yn  and  y will  be  based  on  (6.3) 
and  the  inverse  Radon  transform  (4.11).  To  see  how  (6.3)  can  be  used 
let  SI  C  S2  and  suppose  that  mappings 


(6.4) 


e  :  n  •+  s2,  n  -*•  0(n) , 
e0  :  ft  -*■  s2,  n  -*■  e0(n) , 


can  be  found  such  that 


(6.5)  0(n)  *  0fl(n)  *  U  *  const,  for  n  e  0, 

and 


(6.6)  0(n)  -  0o(n)  =  cy  n  for  n  £  ft 
where 

(6.7)  cu  =  |0(n)  -  80(n)|  -  /2(1-  y)  for  n  e  si  . 

Making  the  substitutions  0  -*•  0(H),  0O  0O  (n)  and  x  -*>  cyT  in  (6.3)  gives 

c3 

(6.8)  h  (x,n)  »  —  e»(c  T,0(p),8o(n))  for  n  6  ft  , 

M  O  M 


and  hence 


(6.9) 


h  (x*n,n)  dn  -  —  e.(c  x  •  n,9(n) ,B0(n))  dn  . 
a  M  2ti  'ft  y 


In  particular,  if  ft  *  S2  then  (4.11)  implies  that  (6.9)  provides  an 
explicit  construction  of  Y^(x)  from  the  scattering  data.  However, 
mappings  9(n)  ,  90(n)  may  only  exist  on  proper  subdomains  ft  c  S2.  In 
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any  case,  if  a  decomposition 


(6.10) 


s2  =  n,  u  n2  u  •••  u  ft 


and  corresponding  mappings 


(6.11) 


9J  :  (L  ->■  S2  . 


S;!  :  flj  -  S2  , 


satisfying  (6.5)  -  (6.7)  can  be  found,  then  one  has 


V*' 


h  (x  •  n,n)  dri 


.2  u 


(6.12) 


!  ' 
3-1 


Q. 


h,,(x  *  n.n)  dn 


c 3  N  , 

=  —  l  e.(c  x  •  n,9:i(n)  ,0-’(n))  dn 

2 it  j=i  Jr  0  u 


If  this  can  be  done  for  two  distinct  values  of  y,  say  Uj  and  u2,  then 
and  yp  can  be  computed  and  yn>  yp  can  be  found  from  the 
equations 


(6.13) 


Y„  +  (U i  -  1) Yp  , 


'V!  n 


Yy2  =  Yn  +  (u*  "  1)Yp  ’ 


The  section  is  concluded  with  a  description  of  a  particular  method  for 
carrying  out  this  program. 

Back  Scattering.  Back  scattering  is  characterized  by  the 


condition  8  *  -80 .  This  may  be  realized  in  the  context  of  equations 
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(6.4)  -  (6.7)  by  defining  Q  -  S2,  and 

6(n)  =  n 

(6.14)  for  all  n  €  S2  . 

eo(n)  =  -n 

Then  (6.5)  -  (6.7)  hold  with  u  =  -1  and  cp  =  c_j  =  2.  Hence  (6.9)  with 
=»  Sz  and  (4.11)  give  the  relation 


(6.15) 


4  f 

Y_,(x)  =  -  J  2  e5(2x  *  n,n,-n)  dn 


In  the  special  case  of  negligible  density  variations,  Yp(x)  =  0.  one  has 
Y_j  =  Yn  and  (6.15)  is  the  one  parameter  imaging  formula  (1.9)  of  the 
introduction. 

Orthogonal  Scattering.  Orthogonal  scattering  is  characterized 
by  the  condition  0  *  0O  =0;  i.e.,  u  =  0.  To  obtain  fields  0(n)  ,  90(n) 
that  satisfy  (6.4)  -  (6.7)  with  U  =  0,  fix  a  geographical  coordinate 
system  in  S2  with  colatitude  a  and  longitude  B  and  define 

(6.16)  0  =  (cos  B  sin  a,  sin  B  sin  a,  cos  a) 

and 
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0O  *  ■  (cos  3  cos  a,  sin  6  cos  a,  -sin  a) 

(6.17) 

-  (cos  B  sin  (a  +  ~) ,  sin  8  sin  (a  +  j)  ,  cos  (a  +  j))  . 

It  is  clear  from  the  identity  6  •  0  ■  1  that  0  •  0O  =0.  Moreover,  a 
short  calculation  gives,  see  Figure  2: 
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(6.18)  ti- —  (9-0o)  *  (cos  6  sin  (a-|),  sin  6  sin  (a-|),  cos  (a-^r)) 


Figure  2.  Graphical  definition  of  0,  0O  and  n- 

For  the  remainder  of  the  discussion  it  will  be  convenient  to  make  the 
change  of  parameter  a  a  +  u/4  and  define 

n  *  (cos  6  Bin  a,  sin  8  sin  a,  cos  a), 

(6.19) ;  0(n)  *  (cos  6  sin  (a+^-)»  sin  8  sin  (a+^),  cos  (a  +  ^)), 

90(n)  •  (cos  8  sin  (o  +  ^p),  sin  8  sin  (a  +  -^),  cos  (a+^r)), 

where 

(6.20)  0  <  a  <  ir  and  0  <  8  <  2ir. 


With  these  choices  n  ranges  over  the  entire  unit  sphere;  l.e.,  ft  •  S*  in 
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(6.4).  The  coordinates  a,  0  are  regular  on  S2  except  at  the  north  and 
south  poles:  a  ■  0  and  a  *  tt.  It  follows  that  9(n)  and  90(n)  are  also 
regular  except  at  these  points  where  they  have  discontinuities.  In 
principle,  these  discontinuities  cause  no  problem  in  evaluating  the 
inverse  Radon  transform.  Hence,  (4.11)  and  (6.9)  with  Q  *  S2  give 


(6.21) 


Y0(x) 


/2 

T T 


Js2  es(^x 


n,8(n) ,0o(n))  dn  , 


where  9(n)  and  90  (n)  are  defined  by  (6.19).  Of  course,  if  one  wishes  to 
avoid  discontinuous  fields  then  the  decomposition  (6.10)  -  t6.12)  may  be 
used  with  different  geographical  coordinate  systems  for  each  component 

Equations  (6.15)  and  (6.21)  define  an  imaging  method  for  the  two 
parameter  problem  because 


(6.22) 


Yp  *  Yo  -  Y_ j »  and 
Yn  '  2Yo  -  Y_j  • 


Of  course,  in  practice  the  integrals  in  (6.15)  and  (6.21)  will  be 
approximated  by  means  of  a  numerical  quadrature  method.  If  the 
quadrature  formula  is 


(6.23) 


M 

F(n)  dn  -  l  F(nk)  Ank 
s2  k-i  k  k 


where  Hj.  n2»*'*»  6  S2  and  Ank  are  suitable  weights  then  the  algorithm 

for  computing  y  and  y0  is 


Y_! (x) 


•«<2x*  wv  Ank 


(6.24) 
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